//rcclust.c
//Programmer: Rory Carmichael
//Date: 1/21/2011
//Purpose: C implementation of pvaluelooper.pl

#include <stdio.h>
#include <unistd.h>
#include <ctype.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <math.h>

float* calc_peak_frequencies(char *filename, int *numcols) {
	int colsize = 0;
	int i = 0;
	int j = 0;
	int first = 1;
	float *freqs;
	int *lengths;
	int *peaks;
	FILE *infile;
	char line[ 10000 ];
	char line2[ 10000 ];
	char delims[] = "\t\n";
	char *result = NULL;
	colsize = *numcols;
	infile = fopen(filename, "rt");
	if ( infile != NULL ) {
		while ( fgets(line, sizeof(line), infile) != NULL) {
			if(first && line != NULL) {
				fgets(line, sizeof(line), infile);
				fgets(line, sizeof(line), infile);
				first = 0;
				memcpy(line2,line,sizeof(char)*10000);
				result = strtok(line, delims);
				while ( result != NULL && (strcmp(result, "\n") != 0)) {
					colsize++;
					result = strtok(NULL, delims);
				}
				lengths = (int *) malloc(colsize*sizeof(int));
				peaks = (int *) malloc(colsize*sizeof(int));
				freqs = (float *) malloc(colsize*sizeof(float));
				for(i=0; i<colsize; i++) {
					lengths[i] = 0;
					peaks[i] = 0;
					freqs[i] = 0.0;
				}
				result = strtok(line2, delims);
				i = 0;
				while ( result != NULL && (strcmp(result, "\n") != 0)) {
					if (strcmp(result, "-") != 0) {
						lengths[i]++;
						if (strcmp(result, "0") == 0) {
							peaks[i]++;
						}
					}
					i++;
					result = strtok(NULL, delims);
				}
			} else {
				result = strtok(line, delims);
				i = 0;
				while ( result != NULL && (strcmp(result, "\n") != 0)) {
					if (strcmp(result, "-") != 0) {
						lengths[i]++;
						if (strcmp(result, "0") == 0) {
							peaks[i]++;
						}
					}
					i++;
					result = strtok(NULL, delims);
				}
			}
		}
		fclose(infile);
		for(i=0;i<colsize;i++) { 
			freqs[i] = ((float) peaks[i]) / ((float) lengths[i]);
		}
		*numcols = colsize;
		return freqs;

	} else {
		printf("Error while opening %s\n", filename);
		perror(filename);
	}
	return NULL;
}

void set_i(int *i) {
	*i = 10;
}

int count_set_bits( uint32_t n ) {
	n = n - (( n >> 1 ) & 0x55555555);
	n = ( n & 0x33333333) + (( n >> 2 ) & 0x33333333);
	return ((n + (n >> 4) & 0xF0F0F0F) * 0x1010101) >> 24;
}

int is_set(int position, uint32_t n) {
	return n & 1<<position;
}

int compare_float(const void *a, const void *b) {
	float fa = *(float *) a;
	float fb = *(float *) b;
	return (int) (fa - fb);
}

int compare_int(const void *a, const void *b) {
	int ia = *(int *) a;
	int ib = *(int *) b;
	return (int) (ia - ib);
}

float *calc_cell_pvalue(int dist, float *pvallist, int numcols) {
	int i = 0;
	float *pcalc = NULL;
	pcalc = (float *) malloc(numcols*sizeof(float));
	for(i=0;i<numcols;i++) {
		pcalc[i] = 1 - pow( (1-pvallist[i]) , (2*dist)+1 );
	}
	return pcalc;
}

char *as_binary(uint32_t num) {
	char *bitstring;
	bitstring = (char *) malloc(33*sizeof(char));
	bitstring[32] = '\0';
	int i = 0;
	while(i < 32) {
		if (num & 1<<i) {
			bitstring[31-i] = '1';
		} else {
			bitstring[31-i] = '0';
		}
		i++;
	}
	return bitstring;
}

float bernoulli_exact(int tgt_num, float *pcalc, int pcalc_len) {
	int sum = 0;
	float microprob = 1.0;
	float semiprob = 0.0;	
	float bitprob;
	uint32_t max = (uint32_t) pow(2, pcalc_len);
	uint32_t i;
	int j;
	for(i=0; i<max; i++) {
		sum = 0;
		microprob = 1.0;
		sum = __builtin_popcount(i);
		for(j=0; j < pcalc_len; j++) {
			if(is_set(j, i)) {
				bitprob = pcalc[j];
			} else {
				bitprob = 1 - pcalc[j];
			}
			microprob *= bitprob;
		}
		if(sum >= tgt_num) {
			semiprob += microprob;
		}
	}
	if(semiprob == 0) { semiprob = 1; }
	return semiprob;
}

float calc_row_pvalue(int numcols, int *distlist, float *pvallist) {
	int i = 0;
	int j = 0;
	float masterprob = 1.0;
	float semiprob = 1.0;
	float *cellprobs;
	qsort(distlist,numcols,sizeof(int),compare_int); 
	for(i = 0; i < numcols; i++) {
		cellprobs = calc_cell_pvalue(distlist[i],pvallist,numcols);
		semiprob = bernoulli_exact(i,cellprobs,numcols);
		if(semiprob < masterprob) {
			masterprob = semiprob;
		}
	}
	return masterprob;
}

int main (int argc, char **argv) {
	FILE *infile;
	int i = 0;
	int j = 0;
	int index = 0;
	float *pvallist = NULL;
	char *filename = NULL;
	char *result = NULL;
	char delims[] = "\t\n";
	char line [ 10000 ];
	int *distlist = NULL;
	int *curdists = NULL;
	float *curpvals = NULL;
	int numcols = 0;
	int gapcount = 0;
	float pval = 1.0;
	char pvalbuf [1000];

	if (argc - index != 2) {
		fprintf(stderr, "No gene file specified, or too many arguments\n");
		exit(1);
	}
	index++;
	filename = argv[index];
	pvallist = calc_peak_frequencies(filename, &numcols);

	distlist = (int *) malloc(numcols*sizeof(int));
	if(pvallist != NULL) {
		infile = fopen(filename, "rt");
		if(infile == NULL) {
			fprintf(stderr, "Failed to open %s\n", filename);
			exit(1);
		}
		fgets(line,sizeof(line),infile);
		//print the geneid
		if (line != NULL) {
//			result = strtok(line, delims);
			fprintf(stdout, "%s\n\n", line);
		}
		fgets(line,sizeof(line),infile);
		while (fgets(line,sizeof(line),infile) != NULL ) {
			result = strtok(line,delims);
			i = 0;
			gapcount = 0;
			while (result != NULL && (strcmp(result, "\n") != 0)) {
				if (i > numcols) {
					fprintf(stderr, "MALFORMED INPUT FILE... %d (with value %s)is too many columns\n", i, result);
					exit(1);
				}
				if (strcmp(result, "-") != 0) {
					sscanf(result, "%d", &distlist[i]);
				} else {
					distlist[i] = -1;
					gapcount++;
				}
				i++;
				result = strtok(NULL, delims);
			}
			curdists = (int *) malloc((numcols-gapcount)*sizeof(int));
			curpvals = (float *) malloc((numcols-gapcount)*sizeof(float));
			i = 0;
			j = 0;
			while(i < numcols-gapcount) {
				if (distlist[j] != -1) {
					curdists[i] = distlist[j];
					curpvals[i] = pvallist[j];
					i++;
				} 
				j++;
			}
			i = 0;
			pval = calc_row_pvalue(numcols-gapcount, curdists, curpvals);
			fprintf(stdout, "%f\n", pval);
		}
	}
	return 0;
}
